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Abstract. The effects of bubbles on the cerebral bloodflow are difficult to quantify. 
We present a model to calculate how cerebral ischaemia is caused by deformable 
gaseous emboli. The model takes into account realistic adhesion forces, fluid dynamical 
considerations, a realistic parameterisaton for the rate of bubble dissolution and the 
effects of buoyancy. We find that neglecting deformability of bubbles leads to a vast 
overestimation of ischaemia. The inclusion of buoyancy effects reduces the proportion 
of the vasculature that becomes compromised, but increases blockage times, thus 
lowering the risk of transient ischaemia but increasing the potential for focal injury. We 
also investigate the number and size of bubbles in a sudden shower of emboli that leads 
to persistent ischaemia capable of neuronal injury. Finally we investigate mitigation 
techniques such as insufflation of the operative area with CO2 and alterations in arterial 
pressure. 



1. Introduction 



Neurological symptoms arising from clinical sources of embolisation can be expected to 
vary considerably depending on embolus size, composition and prevalence [Chung et al. 



2007|. Large solid emboli (consisting of plaque or thrombus) usually form in the heart 
or large arteries such as the aorta or carotid, and once released into the circulation can 
have devastating consequences, travelling through the blood stream to become lodged 
in smaller arteries supplying the cortex. Embolisation can be asymptomatic, or lead to 
a spectrum of neurological disorders ranging from a subtle decline in neuropsychological 
function, or transient ischaemic events (such as decompression sickness), to pronounced 
focal deficits characteristic of severe stroke. 

During invasive procedures such as open heart surgery, patients typically experience 
thousands of gaseous emboli with various clinical implications (Rodriguez et al., 2005 



Moody et al., 1990 . Bubbles can also be observed in other scenarios as a cause of 



decompression sickness Barak and Katz, 2005 and in patients fitted with mechanical 



heart valves Georgiadis et al., 20051. Depending on the number, size and embolisation 
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rate, embolisation of bubbles clearly produces various sequelae and it is the aim of this 
work to understand where the boundaries between symptoms occur. 

To attempt to understand the outcomes from these varied sources of embolisation 
in a universal way, we previously introduced a model for the forecast of ischaemic stroke 
[Chung et al. , 2007, Hague and Chung, 2009 . In that model, emboli were assumed to be 



rigid, and therefore to block any vessels of similar size. The model also assumed that all 
pressure was dropped over the arterioles and the flow at any bifurcation was determined 
only by the number of arterioles receiving flow downstream of the bifurcation. Our 
stroke model was used to discuss various clinical scenarios of embolisation, and revealed a 



nonlinear relationship between embolus properties and cerebrovascular blockages Chung 



et al., 2007]. Previous clinical studies of embolisation revealed a positive correlation 
between embolus prevalence and adverse clinical outcome during some stages of surgery, 
but in many clinical situations the presence of emboli alone is insuflicient to provide a 



2005] . 



, 2006, 


Brown et al. , 


2000, 


Stroobant et al. 



There are a number of differences between solid emboli (which may be almost 
incompressible) and gaseous emboli which compress easily [Branger and Eckmann , 1999 



that have potential to affect the locations in the vasculature where emboli become 
lodged. Gaseous emboli have a propensity to split and deform as they move through 
the vasculature and only block arteries when the surface area in contact with the walls 



becomes sufficiently large to oppose motion through stiction (Suzuki and Eckmann 
|200 3|), whereas solid emboli block arteries when they encounter vessels of similar size. 
There is another significant difference between solid and gaseous emboli, which is their 
high buoyancy compared to thrombus which is approximately neutrally buoyant in the 
blood. Therefore a number of extensions are needed to properly model the motion of 
gaseous emboli through the vascular tree: (1) the emboli should be deformable (2) an 
approximation for stiction should be included, (3) an iterative fiuid dynamical analysis 
is required to estimate the pressure drop in the whole tree, and (4) in a 3D vasculature 
an estimate of buoyancy effects would be useful. 

In this paper, we discuss an extended model which includes features of gaseous 
emboli such as buoyancy and deformation. We start by introducing an extended stroke 
forecast model (section [2]) and then present detailed results for prediction of transient 
and permanent ischaemic injury caused by showers of gaseous emboli, contrasting the 
results with those for solid emboli, and probing the effects of buoyancy (section [3]) . A 
summary can be found in section |4| 



2. Model 



2.1. Gaseous emboli 



Gaseous emboli differ from solid thrombus and plaque since they are deformable, and 
also dissolve at a faster rate. By fitting to in vivo data from [Branger and Eckmann 
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Figure 1. Schematic of a deformable embolus and forces leading to stiction 



1999 it becomes possible to parameterise the rate of embolus dissolution, yielding the 
following function: 

V{t) = c{a - tf (1) 

where the volume at t = corresponds to an embolus of size V{0) = 2.8 nl. The 
parameters for an air embolus are a = {V{0) /cY^^ minutes (which is the lifetime of a 
2.8 nl embolus), b = 2.00539 ± 0.04699 and c = 0.00463395 ± 0.0008206 nl/min^ CO2 
emboli dissolve around 20 times faster ( Chaudhuri and Marasco| [201 1|) so (representing 



CO2 parameters as primed variables) d = 20^ x = b and a' = {V{0)/c^y^^. 

In order to consider deformability and stiction we assume that a gaseous embolus 
in a vessel of similar size deforms to a sausage-like shape (shown schematically in figure 
[1]) . This shape is a common finding in the arteries of patients following cardiac surgery. 



where autopsy reveals large numbers of sausage-like arteriole dilatations (Moody et al. 
[1990] ). The surface area of the embolus touching the side of the vessel is computed 
by correcting for a domed end with the same radius as the vessel. The length of the 
cylindrical part of the embolus is L = (V^mb ~47ri?^ggggj/3)/7ri?^ggggj, and the surface area 
touching the side is 27ri?vessei^- 

An embolus will stop moving when the pressure drop over a stationary embolus is 
insufficient to overcome stiction, i.e. Ap/ttRI^^^^^ < 2X(V^mb ~ ^^Rlessei/^) / -^yessei- The 



coefficient of stiction has been measured, and is K = lONm ^ (Suzuki and Eckmann 
[2003] ). The pressure drop for the stationary embolus is equal to the difference between 
the pressure at the bifurcation upstream from the embolus and the capillary network 
(since the blood in vessels downstream of the embolus is not moving, and therefore there 
is no additional pressure drop over the downstream vessels). 

The need to assign values to the pressures at bifurcations in the arterial tree is an 
added complication of dealing with deformable gaseous emboli, and means that at least 
a basic fluid dynamical analysis of flows in the tree is needed. This will be discussed 
below. 

2.2. Bifurcating tree 

The assumption of a bifurcating tree is justifled by the analysis of high-resolution images 
of the human cortex (Cassot et al. [2006J ) which show that 94% of branches in the 



cerebral vasculature consist of bifurcations [|j 
J Of the remaining nodes, 4% are trifurcations, 1% simple nodes and 0.5% have 4 or more daughters 
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Figure 2. Schematic of the recursive procedure for computing pressures, flows and 
resistances. 



2.3. Fluid dynamics and recursive computation 

When considering the stiction of gas bubbles, it is essential to know the pressure drop 
across the embolus. Since the bifurcating tree in the model has a very large number of 
nodes, and the pressures need to be computed whenever an embolus moves, it is essential 
that a simplified fiuid dynamics scheme is used so that the computational time is not 
excessive. We treat the fluid flow through the tree as Poiseuille flow. The pressure drop 
and effective resistance can be treated using electrical circuit analogues. Thus, parallel 
resistances can be rewritten as a single resistor, and this can be repeated recursively up 
the tree from the end arterioles to the parent node to compute all flows and pressures 
as an order operation, where is the number of vessels in the tree. This process 
is summarised in flgure [2| and leads to rapid computations which are far faster than 
using matrix inversion or by directly solving simultaneous equations. In the following, 
the root node radius of the tree is 1mm. 

Panel (a) shows the initial bifurcating tree. In panel (b), the smallest vessels have 
been rewritten as a single resistance. In panel (c) the resistances in series have been 
simplifled, and in panel (d) the recursive step has been applied to the tree that resulted 
in (c). Clearly this procedure can be applied to any size of bifurcating tree. It is a 
simple matter to compute all flows and pressures at step (b) as there is effectively a 
potential divider in each branch at that stage. During the recursive procedure, these 
flows and pressures are stored, and the stored pressures and flows can be used for all 
calculations until any new blockages are introduced or existing blockages are freed as 
emboli dissolve. 
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2.4' Buoyancy 

The effects of buoyancy are emulated by introducing preferential directions at 
bifurcations. An additional weighting, wa = (1+^^ cos(0))/2 is included in the weighted 
probability to account for buoyancy. This type of weighting is consistent with the results 
in Eshpuniyani et al. 2005 . Once the additional weighting factors are introduced, the 
probability that an embolus travels in direction A at a bifurcation is. 

Pa = wa/a/ {wa/a + wb/b) (2) 

with Pb = 1 — Pa and where /a {/b) is the flow in the a (b) direction. Ag is a 
parameter which is varied between and 1, with Ag = representing no correction 
due to buoyancy (where t(; = 1/2 as in the previous version of the model) and Ag = 1 
an extreme correction. We note that this form for the weighting is ad-hoc^ but we 
believe that it contains the essence of the corrections that will occur when emboli are 
buoyant. In this work we will use it as a tuning parameter to understand some features 
of buoyancy and future work will measure the precise form of the probabilities for a wide 
range of flow ratios. In this work, is a random parameter assigned to each bifurcation 
at time t = to each instance of the ensemble, but in future work it will be the angle 
established from in-silico tree growth. Incidentally, this form may also mimic more 
realistic bifurcating trees, where large trunks have small branches. 

2.5. Algorithm 

The algorithm begins by calculating flows, pressures and resistances for an empty tree 



(using the procedure in section 2.3). It then proceeds as follows: 

(i) On any time step, an embolus may be created in the root node of the tree with 
probability P^-Ar with size randomly chosen between and Tmax- Here is the 
probability per unit time to create an embolus and At = Is is the length of the 
time step. 

(ii) All emboli dissolve leading to a reduction in radius during each time step according 



to the parameterisation in Sec. |2.1[ Completely dissolved emboli are removed from 
the simulation. If the reduction in radius causes a change in the blockage state of 
the tree, flows and pressures are calculated again. 

(iii) The emboli move according to the following rules: 

(a) If the pressure behind the deformed embolus is insuflicient to overcome the 



stiction force exerted by the embolus (according to the theory in section 2.1) 
it does not move. 

(b) If all arterioles downstream are blocked, the embolus may not move since there 
is no flow. 

(c) If the embolus is smaller than the current node, and there is flow downstream. 

1. The flows in directions A and B are determined by solving the simplifled 
fluid dynamics. 
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2. The embolus then moves in direction A with probabihty Pa = 
fA'^A/ifAUJA + IbUJb)- Otherwise, it moves in direction B. 

3. Results 

Here, we investigate total blockage, and the percentage of arterioles blocked for at least 
5 minutes, 10 minutes and 15 minutes as these are the timescales over which we expect 
neuronal changes to occur. As part of these simulations the size range of incoming 
emboli, buoyancy, the total number of emboli in a shower, and the rate of embolisation 
are varied. We also investigate the difference between deformable gaseous emboli and 
non-deformable solid emboli. Embolisation starts at time t = 10s and finishes after the 
number of emboli to be explored has been delivered. The embolisation is steady, and 
occurs at a rate of 1/s and 0.5/s. In this way, we can determine if a short period of 
rapid embolisation is more damaging than chronic embolisation over a longer period. 
In clinical scenarios, the embolisation rate can be higher than this, but it becomes 
difficult to differentiate between multiple emboli. For embolisation showers which have 
a significantly shorter time duration than the embolus dissolve time, the outcome from 
the showers is expected to be similar. An ensemble of 100 simulations is generated, with 
the weighting coefficients for buoyant emboli chosen randomly in each instance. 

Figure |3] shows the percentage of nodes blocked on each timestep. In the simulations 



shown, emboli have sizes ranging between to 0.5 mm. As in our original model (Chung 



et al. 2007| ), after embolisation begins the proportion of blocked nodes quickly reaches 



a dynamic equilibrium where as many emboli are dissolving and leaving the system as 
entering the vasculature. There is a small difference between the height of the plateau as 
buoyancy effects are turned on in the simulation. Buoyancy effects reduce the proportion 
of blocked nodes because certain paths through the tree become more probable if the 
emboli are buoyant. Once embolisation stops, the emboli rapidly dissolve and blockages 
are cleared. 

There is a major difference between the percentage of blocked end arterioles when 
emboli are assumed to be deformable and rigid. If emboli are assumed to be rigid, the 
proportion of blocked nodes rapidly increasing to 100% for the parameters shown. In 
contrast, approximately 5% of end arterioles become blocked when deformable emboli 
are introduced into the model. This shows that deformability must be taken into account 
in models of gaseous embolisation. 

Knowledge of the percentage of instantaneous blockages is of limited diagnostic 
use, since it does not contain any information about the length of time that a particular 
arteriole has received zero ffow. If ffow is absent for a sufficiently long time, neuronal 
changes begin to occur (<1 min) and surrounding tissue becomes ischaemic (<5 mins). 
Eventually these changes become irreversible leading to necrosis and cerebral infarcts. 
Therefore, we also measure the proportion of nodes that have been ischaemic for at least 
5 minutes, 10 minutes etc. 

An estimate for the number of nodes that have been ischaemic for at least 10 minutes 
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Figure 3. Percentage of nodes blocked vs time giving an indication of transient 
effects. After embolisation begins the proportion of blocked nodes quickly reaches a 
dynamic equilibrium where as many emboli are dissolving and leaving the system as 
entering the vasculature. There is a small difference between the height of the plateau 
as gravitational effects are turned on in the simulation, and gravitational effects also 
reduce the proportion of blocked nodes because certain paths through the tree become 
more probable if the emboli are buoyant. Once embolisation stops, the emboli rapidly 
dissolve and total blockage returns to zero. N.B. The graph is plotted on a logarithmic 
scale. 
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Figure 4. Percentage of individual nodes compromised for at least 10 minutes vs time. 
Again several values for the coefficient Ag are shown, and the results for non-deformable 
emboli are included for comparison. The proportion of nodes that have been ischaemic 
for at least 10 minutes at some time during the simulation increases monotonically 
with time and the number of emboli in the shower. It should be noticed that these 
proportions do not increase linearly, presumably because the same nodes can become 
ischaemic at a later time in the simulation. This non-linearity is particularly clear 
for Ag = 1 since emboli are weighted in the same direction. The effect of buoyancy 
is to increase the proportion of ischaemic nodes, which is also a consequence of this 
weighting. Again, a lack of deformation in emboli leads to near 100% ischaemia. N.B. 
Results are plotted on a logarithmic scale. 
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Figure 5. The proportion of end arterioles receiving no supply for 5 minutes when (a) 
buoyancy is turned off and 1 embolus is produced per second (b) buoyancy is turned 
off and 0.5 an embolus is produced per second (c) Ag = 0.5 and (d) Ag = 1. The 
shape of the graph is similar in all cases, but the proportion of blocked nodes varies 
dramatically (the colour scales vary between panels). Between panels (a) and (b), 
the only difference is the rate of embolisation, showing that it is the total number of 
emboli rather than the exact pattern of embolisation that is most important when the 
timescale of embolisation is lower than the dissolve time. Panels (c) and (d) show 
that buoyancy reduces the maximum blockage for the parameters sampled to 25% for 
Ag = 0.5 and to 10% for Ag = 1. 

is plotted vs time in figure [4[ Again several values for the coefficient Ag are shown, and 
the results for non-deformable emboli are included for comparison. The proportion 
of nodes that have been ischaemic for at least 10 minutes at some time during the 
simulation increases monotonically with time and the number of emboli in the shower. 
It should be noticed that these proportions do not increase linearly, presumably because 
the same nodes can become ischaemic at a later time in the simulation. This non- 
linearity is particularly clear for Ag = 1 since emboli are weighted in the same direction. 
The effect of buoyancy is to increase the proportion of nodes that are ischaemic for 
a long time. A lack of deformation in emboli leads to near 100% ischaemia. This 
overestimation can also be seen in Fig. [6| 

Figure [5] shows the percentage of arterioles subjected to at least 5 minutes of 
ischaemia for gaseous emboli. Increasing the effects of buoyancy decreases the total 
percentage of arterioles blocked as the paths of emboli through the tree are restricted 
by the weighting. Overall, an increase in the maximum size r^ax of emboli increases the 
proportion of arterioles affected by ischaemia. Likewise, increasing the number of emboli 
in the shower is predicted to increase the proportion of the tree affected by ischaemia. If 
no blockage is registered, then the cell is black. These calculations demonstrate that a 
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Figure 6. Comparison between tiie proportion of nodes receiving no blood supply 
for 5 minutes when (a) emboli are non-deformable and (b) emboli are deformable. 
Buoyancy effects are not considered. Neglect of embolus deformation in the model 
leads to a huge overestimation of the proportion of blocked nodes. 
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Figure 7. The proportion of nodes receiving no supply for 10 minutes when (a) 
buoyancy is turned off and 1 embolus is produced per second (b) buoyancy is turned 
off and 0.5 an embolus is produced per second (c) Ag = 0.5 and (d) Ag = 1. 



smaller number or size of emboli is required to cause ischaemia when emboli are buoyant. 

To highlight the differences between deformable and non-deformable emboli, figure 
[g] shows the proportion of nodes blocked for more than 5 minutes when emboli are (a) 
rigid and (b) deformable. There are significantly more blockages in the case of non- 
deformable emboli (in both cases, the same clearance time is used), and this difference 
persists over all timescales for the ischaemia. 

It is of interest to see how ischaemia is affected over longer timescales. Figure [7| 
(figurejs]) shows the percentage of arterioles experiencing at least 10 minutes (15 minutes) 
of ischaemia. The proportion of arterioles blocked for a longer period of time is around 
one order of magnitude smaller than the proportion blocked for at least 5 minutes, since 
the clearance time for a 0.5 mm gaseous embolus is just under 10 minutes [Branger and 
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Figure 8. The proportion of nodes receiving no supply for 15 minutes when (a) 
buoyancy is turned off and 1 embolus is produced per second (b) buoyancy is turned 
off and 0.5 an embolus is produced per second {c) Ag = 0.5 and (d) Ag = 1. 



Eckmann, 1999] and obstruction from a single embolus drops exponentially with size. 



Of particular interest is the observation that buoyant emboli produce more persistent 
blockages (up to 1.8% of arterioles blocked for at least 10 minutes for 1500 gaseous 
emboli and 0.6% for at least 15 minutes). This is a reversal of the short timescale 
reduction in blockage, and demonstrates that if emboli have preferred paths through 
the tree (caused either by topography or buoyancy) their effects are more likely to 
translate to tissue ischaemia and neuronal injury. 

Finally, we discuss two possible mitigation techniques for dealing with gaseous 
emboli. The first is increase of blood pressure, and the second is the use of CO2 during 
surgery. Results for these scenarios can be found in Fig. [9] for emboli of radius up to 
0.5mm, and pressures of 100 mmHg and 50 mmHg. Inclusion of pure CO2 dramatically 
reduces the blockage, with less than 1% of arterioles blocked on any time step (panel 
(a)), and with no nodes experiencing ischaemia for > 5mins, which demonstrates that 
under ideal conditions, CO2 insufflation should be effective in reducing the contribution 
of air bubbles to neurological impairment during cardiac surgery. Another interesting 
result is the effect of decreased pressure. A drop of the input pressure by a factor of 
2 leads to a near doubling in both the number of blocked nodes during embolisation, 
and in the number of nodes suffering ischaemia for more than 5 minutes. While there 
may be other good clinical reasons why a lower pressure might be useful on bypass, the 
results here indicate that maintaining higher pressures may reduce stroke symptoms. 
The origin of the drop is that the higher pressure can push gaseous emboli further into 
the tree before they stick, reducing the blockage since the number of nodes downstream 
from a blockage halves at each bifurcation. 
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Figure 9. Effect of halving pressure or adding CO2 on (a) the percentage of blocked 
nodes and (b) the percentage of nodes that have been ischaemic for at least 5 minutes 
during the simulation. The results demonstrate that, under ideal conditions, CO2 
insufflation should be effective in reducing the contribution of gaseous emboli to 
neurological injury during cardiac surgery. A drop of the input pressure by a factor of 
2 leads to a near doubling in both the number of blocked nodes during embolisation, 
and in the number of nodes suffering ischaemia for longer than 5 minutes. N.B. for 
the CO2 bubbles, ischaemia did not exceed 5 minutes. 



4. Summary 

We have investigated ischaemia due to gaseous emboh using an extended model. Our 
model includes the deformability of gaseous emboli, fluid dynamical considerations and 
a basic description of embolus buoyancy. We have examined the effects of buoyancy, 
embolus composition, embolus rate, embolus size, input pressure and the number of 
emboli in a shower on ischaemia. 

We found that deformations of gas bubbles must be considered or blockages in our 
simulations are vastly overestimated. We also found that gravitational effects reduced 
the proportion of blocked nodes at a speciflc time, but slightly increased the proportion 
of nodes that were ischaemic for an extended period of time. Finally, we discussed two 
possible mitigation techniques for dealing with gaseous emboli, confirming that under 
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ideal conditions, CO2 insufflation should be effective in reducing strokes. We also found 
that maintaining higher blood pressures in patients on bypass may reduce cognitive 
decline, with the caveat that it may lead to other clinical complications. 

We note that the model still has some limitations. Pulsatile flow has not been 
considered so far, but flow in the microvasculature is in a near steady state, so the 
effects of pulsatile flow are likely to be small. We expect to include windkessel equations 
in a future model to handle pulsatile flow. In our model, gaseous emboli do not split 
at bifurcations. However, bubble splitting is only likely to affect bifurcations angled 
within, or close to, the horizontal plane. Finally, our tree is completely symmetric. We 
are currently in the process of growing realistic cerebral vasculatures in-silico. 
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